Raghav K. Chhetri

07-21-2021

• Discrete Fourier Transform (DFT) is a way to transform a discrete signal from time to frequency domain by summing over a sequence of sine and cosine waves. Fast Fourier Transform (FFT) is a series of algorithms that compute DFT by factorizing the matrix into a product of sparse matrices, thereby reducing the complexity and speed of computation enormously

There are multiple ways to compute FFT in Python:

1. numpy.fft
2. scipy.fftpack -- Considered legacy. SciPy recommends using scipy.fft
3. scipy.fft
4. pyfftw
5. cupy.fft -- NVIDIA's library that evaluates FFTs on GPUs
6. torch.FFT

• 2D FFT speed test shows torch.fft as the fastest option
• scipy.fft is slightly faster than numpy.fft, so using that in this notebook

import sys
print('Environment: {}'.format(sys.exec_prefix))
print('Executable: {}'.format(sys.executable))
print('Version: {}'.format(sys.version))

Environment: D:\Anaconda3\envs\denoising
Executable: D:\Anaconda3\envs\denoising\python.exe
Version: 3.7.10 | packaged by conda-forge | (default, Feb 19 2021, 15:37:01) [MSC v.1916 64 bit (AMD64)]

import numpy as np

import scipy.misc
from scipy import ndimage, signal
from scipy.fft import fft, ifft, fftn, ifftn, fftshift, ifftshift, fftfreq, rfft, irfft, rfftfreq

import skimage
from skimage.util import img_as_float
from skimage.data import camera, gravel, checkerboard
from skimage.filters import difference_of_gaussians, window

import matplotlib
from matplotlib import cm
import matplotlib.pyplot as plt
#plt.style.use('seaborn-poster')
%matplotlib inline

print('numpy: {}'.format(np.__version__))
print('scipy: {}'.format(scipy.__version__))
print('skimage: {}'.format(skimage.__version__))
print('matplotlib: {}'.format(matplotlib.__version__))

numpy: 1.18.5
scipy: 1.6.0
skimage: 0.17.2
matplotlib: 3.3.4

---

1D FFT (for time-series data)

General comments:

fft takes longer if npnts is not a multiple of 2 and even longer if it is a prime

The frequency axis has been defined in different ways in fft_powerspectrum_plot below:

• Method 1 (Frequency mirroring above Nyquist):\ Manually define the frequency axis in this range, [0, Sampling Rate]. The FFT Amplitude of the signal now has frequency peaks that are mirrored in the [Nyquist, Sampling Rate] too. Ignore the mirrored peaks and only consider this range, [0 to Nyquist frequency]. Recall that Nyquist frequency is half of the Sampling Rate.

• Method 2 (Frequency mirroring about zero):\ Define the frequency axis using fftfreq which centers the output of fft(). The frequency axis now extends from [-Nyquist, Nyquist] and the peaks are mirrored in the negative axis. Again, ignore the negative frequency and only consider this range, [0 to Nyquist frequency]

• Method 3 (No mirroring-- only positive frequencies):\ The frequency spectrum that fft() outputted was reflected about the y-axis, which is caused by inputting real numbers instead of complex numbers to fft(). This symmetry can be taken advantage of to make Fourier transform faster by computing only the positive frequencies, which is what rfft() does.

def fft_powerspectrum_plot(x, t, dt, npnts, freq_show, method=2, denoise=False, cutoff=0):   
    ''' 
    - plot original signal
    - compute FFT and power spectrum (power per frequency)
    - plot power spectrum
    - plot recovered signal after fft-> ifft
     
    x:         1D signal
    t:         time vector
    dt:        sampling interval 
    npnts:     number of time points
    freq_show: upper range of frequency to plot
    method:    1/2/3
    denoise:   filter out noise
    cutoff:    denoise below the cutoff amplitude 
    
    Call: plot_signal_fftamplitude(x, t, dt, npnts, 200, method=1, denoise=True, cutoff=3)
    '''
    if method == 1:
        X = fft(x)
        n = np.arange(npnts)
        T = npnts*dt
        freq = n/T 
        title = 'Frequency mirroring above Nyquist'
    elif method == 2:
        X = fft(x)
        freq = fftfreq(npnts, dt)
        title = 'Frequency mirroring about zero'
    elif method == 3:   
        X = rfft(x)
        freq = rfftfreq(npnts, dt)
        title = 'No mirroring: only positive frequencies'
    
    powerspect = 2*np.abs(X)/npnts   
    # Note: Normalized as 2*np.abs(X)/npnts instead of simply np.abs(X)
    # returns the actual amplitude values of the sine and cosine functions
    # instead of some arbitrarily-scaled values
    
    if denoise:
        powerspect = powerspect * (powerspect > cutoff) # Zero all frequencies with small power
        X = X * (powerspect > cutoff) # Zero small Fourier coefficients        
        #To further zero out a peak at zero frequency -- occurs if noise is 'rand' instead of 'randn'
        #powerspect = powerspect * (powerspect < 10) 
        #X = X * (powerspect < 10)

    # PLOT
    plt.figure(figsize = (9, 7))
    plt.subplot(311)
    plt.plot(t, x, 'k',label='original')
    plt.xlabel('Time (s)')
    plt.ylabel('Amplitude')
    plt.legend()
    
    plt.subplot(312)
    plt.title(title)
    #plt.stem(freq, np.abs(X),'c', markerfmt=" ", basefmt="-b")
    plt.stem(freq, powerspect,'c', markerfmt=" ", basefmt="-b")
    plt.xlabel('Freq (Hz)')
    plt.ylabel('FFT Amplitude')
    if method == 2:
        plt.xlim(-freq_show,freq_show)
    else:
        plt.xlim(0,freq_show)

    plt.subplot(313)
    if method == 3:
        plt.plot(t, irfft(X), 'k--',label='recovered')
    else:
        plt.plot(t, ifft(X), 'k--',label='recovered')
    plt.xlabel('Time (s)')
    plt.ylabel('Amplitude')
    plt.legend()

    plt.tight_layout()
    plt.show()
    
    print('Number of time points:', npnts, 'points')
    print('Number of points in frequency range:', len(freq), 'points')
    print('Frequency range:', min(freq), max(freq), 'Hz')

Let's generate a 1D signal and add some noise to it:

srate = 1000 # Sampling rate in Hz
dt = 1/srate
duration = 1 # in seconds
t = np.arange(0,duration,dt) # Time Vector in seconds
npnts  = len(t)

# 1D signal
freq1 = 10.0
x = 2*np.sin(2*np.pi*freq1*t)

freq2 = 70.0
x += 4*np.cos(2*np.pi*freq2*t)

freq3 = 170.0  
x += 6*np.cos(2*np.pi*freq3*t+np.pi/4)

noise = 6.0*np.random.randn(npnts)
x += noise

Let's plot to see what the signal looks like and what its power spectrum looks like (no denoising yet):

fft_powerspectrum_plot(x, t, dt, npnts, 200, 1)

D:\Anaconda3\envs\denoising\lib\site-packages\numpy\core\_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Number of time points: 1000 points
Number of points in frequency range: 1000 points
Frequency range: 0.0 999.0 Hz

Now, let's denoise frequencies with power spectrum below 1:

fft_powerspectrum_plot(x, t, dt, npnts, 200, 1, True, 1)

D:\Anaconda3\envs\denoising\lib\site-packages\numpy\core\_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Number of time points: 1000 points
Number of points in frequency range: 1000 points
Frequency range: 0.0 999.0 Hz

If we plot the frequency axis to show mirrored values in the negative axis too, this is what it looks like:

fft_powerspectrum_plot(x, t, dt, npnts, 200, 2, True, 1)

D:\Anaconda3\envs\denoising\lib\site-packages\numpy\core\_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Number of time points: 1000 points
Number of points in frequency range: 1000 points
Frequency range: -500.0 499.0 Hz

Or, using method 3, which uses rfft instead of fft:

fft_powerspectrum_plot(x, t, dt, npnts, 200, 3, True, 1)

Number of time points: 1000 points
Number of points in frequency range: 501 points
Frequency range: 0.0 500.0 Hz

---

Spurious frequency peaks in Amplitude Spectrum

Let's change the frequencies of the signal to non-integer values. When we do so, we see that the FFT amplitude shows up with spurious frequency peaks

srate = 1000 # Sampling rate in Hz
dt = 1/srate
duration = 1 # in seconds
t = np.arange(0,duration,dt) # Time Vector in seconds
npnts  = len(t)

# 1D signal
freq1 = 10.1
x = 2*np.sin(2*np.pi*freq1*t)

freq2 = 70.5
x += 4*np.cos(2*np.pi*freq2*t)

freq3 = 170.7  
x += 6* np.cos(2*np.pi*freq3*t+np.pi/4)

noise = 6.0*np.random.randn(npnts)
x += noise

fft_powerspectrum_plot(x, t, dt, npnts, 200)

D:\Anaconda3\envs\denoising\lib\site-packages\numpy\core\_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Number of time points: 1000 points
Number of points in frequency range: 1000 points
Frequency range: -500.0 499.0 Hz

One way to mitigate these additional frequency components is to use a longer signal, as shown below:

srate = 1000 # Sampling rate in Hz
dt = 1/srate
duration = 10 # in seconds
t = np.arange(0,duration,dt) # Time Vector in seconds
npnts  = len(t)

# 1D signal
freq1 = 10.1
x = 2*np.sin(2*np.pi*freq1*t)

freq2 = 70.5
x += 4*np.cos(2*np.pi*freq2*t)

freq3 = 170.7  
x += 6*np.cos(2*np.pi*freq3*t+np.pi/4)

noise = 6.0*np.random.randn(npnts)
x += noise

fft_powerspectrum_plot(x, t, dt, npnts, 200)

D:\Anaconda3\envs\denoising\lib\site-packages\numpy\core\_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Number of time points: 10000 points
Number of points in frequency range: 10000 points
Frequency range: -500.0 499.90000000000003 Hz

Another option is to denoise lower frequencies like we did earlier

---

2D FFT (working with images)

1. Blurring an image with 2D FFT

Adapted from the following numpy example

tmp = scipy.misc.face(gray=True) #Racoon face from SciPy
image = tmp[128:640, 256:768]
image.shape

(512, 512)

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('Original image')
plt.axis('off')
plt.imshow(image,cmap='gray');

# Take the 2-dimensional FFT and center the frequencies
ftimage = fftn(image)
ftimage = fftshift(ftimage)

plt.subplot(122)
#plt.imshow(np.abs(ftimage))
plt.imshow(np.log(np.abs(ftimage)), cmap='gray')
plt.title('FT(image)')
plt.axis('off')
plt.show()

# Build and apply a Gaussian filter
sigmax, sigmay = 7, 13
size = image.shape[0]
cy, cx = size/2, size/2
x = np.linspace(0, size, size)
y = np.linspace(0, size, size)
X, Y = np.meshgrid(x, y)
gmask = np.exp(-(((X-cx)/sigmax)**2 + ((Y-cy)/sigmay)**2))
#Note: mask is in Fourier space; same size as FT(image)

plt.figure(figsize = (12, 6))
plt.subplot(131)
plt.imshow(gmask, cmap='gray')
plt.title('Gaussian mask')

ftimagep = ftimage * gmask
plt.subplot(132)
plt.imshow(np.log(np.abs(ftimagep)),cmap='gray')
plt.title('FT(image)*mask')
plt.axis('off')

# Finally, take the inverse transform and show the blurred image
imagep = np.fft.ifftn(ftimagep)

plt.subplot(133)
plt.imshow(np.abs(imagep), cmap='gray')
plt.title('Blurred image')
plt.axis('off')
plt.show()

D:\Anaconda3\envs\denoising\lib\site-packages\ipykernel_launcher.py:35: RuntimeWarning: divide by zero encountered in log

---

2. Gaussian blur implemented using FFT convolution: fftconvolve

#From https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.fftconvolve.html
window_x = signal.windows.gaussian(101, std=7)
window_y = signal.windows.gaussian(101, std=13)

plt.plot(window_x, label= "$\sigma_x$=7")
plt.plot(window_y, label= "$\sigma_y$=13")
plt.legend(loc="upper right")
plt.title("Gaussian window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")

kernel = np.outer(window_x, window_y)
#Note: Kernel is in the same coordinate as the image; Its size is smaller compared to image

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.imshow(kernel, cmap='gray')
plt.title('Gaussian Kernel')

blurred = signal.fftconvolve(image, kernel, mode='same')

plt.subplot(122)
plt.imshow(blurred, cmap='gray')
plt.title('Blurred image')
plt.axis('off')
plt.show()

Notice the dark borders around the image, due to zero-padding beyond its boundaries. scipy.ndimage.gaussian_filter gets rid of this artifact.

---

3. Gaussian blur using SciPy.ndimage

Adapted from the following SciPy lecture

# Image Denoising: Gaussian filter, Median filter
noisy = image + 0.4*image.std()*np.random.random(image.shape)

gauss_denoised = ndimage.gaussian_filter(noisy, sigma = 11)
med_denoised = ndimage.median_filter(noisy, size = 9)

plt.figure(figsize=(12,4))

plt.subplot(131)
plt.imshow(noisy, cmap=plt.cm.gray, vmin=40, vmax=220)
plt.axis('off')
plt.title('noisy', fontsize=20)
plt.subplot(132)
plt.imshow(gauss_denoised, cmap=plt.cm.gray, vmin=40, vmax=220)
plt.axis('off')
plt.title('Gaussian filter', fontsize=20)
plt.subplot(133)
plt.imshow(med_denoised, cmap=plt.cm.gray, vmin=40, vmax=220)
plt.axis('off')
plt.title('Median filter', fontsize=20)

plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, bottom=0, left=0, right=1)
plt.show()

Note that a Gaussian filter smooths out the noise and the edges. Median filter averages the noise and preserves the edges better

---

Band-pass filtering by Difference of Gaussians

Adapted from the following SciKit-Image documentation

Band-pass filters attenuate signal frequencies outside of a range (band) of interest. In image analysis, they can be used to denoise images while at the same time reducing low-frequency artifacts such a uneven illumination. Band-pass filters can be used to find image features such as blobs and edges.

One method for applying band-pass filters to images is to subtract an image blurred with a Gaussian kernel from a less-blurred image. This example shows two applications of the Difference of Gaussians approach for band-pass filtering.

More on DoG here

1. Denoise image and reduce shadows

image = gravel()
wimage = image * window('hann', image.shape)  # window image to improve FFT
filtered_image = difference_of_gaussians(image, 1, 12)
filtered_wimage = filtered_image * window('hann', image.shape)
im_f_mag = fftshift(np.abs(fftn(wimage)))
fim_f_mag = fftshift(np.abs(fftn(filtered_wimage)))

fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(11, 12))
ax[0, 0].imshow(image, cmap='gray')
ax[0, 0].set_title('Original Image')
ax[0, 0].axis('off')
ax[0, 1].imshow(np.log(im_f_mag), cmap='magma')
ax[0, 1].set_title('FFT Magnitude (log)')
ax[0, 1].axis('off')
ax[1, 0].imshow(filtered_image, cmap='gray')
ax[1, 0].set_title('Filtered Image')
ax[1, 0].axis('off')
ax[1, 1].imshow(np.log(fim_f_mag), cmap='magma')
ax[1, 1].set_title('FFT Magnitude (log)')
ax[1, 1].axis('off')
plt.show()

2. Enhance edges in an image

image = camera()
#image = noisy
wimage = image * window('hann', image.shape)  # window image to improve FFT
filtered_image = difference_of_gaussians(image, 1.5)
filtered_wimage = filtered_image * window('hann', image.shape)
im_f_mag = fftshift(np.abs(fftn(wimage)))
fim_f_mag = fftshift(np.abs(fftn(filtered_wimage)))

fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(11, 12))
ax[0, 0].imshow(image, cmap='gray')
ax[0, 0].set_title('Original Image')
ax[0, 0].axis('off')
ax[0, 1].imshow(np.log(im_f_mag), cmap='magma')
ax[0, 1].set_title('Original FFT Magnitude (log)')
ax[0, 1].axis('off')
ax[1, 0].imshow(filtered_image, cmap='gray')
ax[1, 0].set_title('Filtered Image')
ax[1, 0].axis('off')
ax[1, 1].imshow(np.log(fim_f_mag), cmap='magma')
ax[1, 1].set_title('Filtered FFT Magnitude (log)')
ax[1, 1].axis('off')
plt.show()

---

Diffraction-limited imaging simulation

• Adopted from Computational Fourier Optics: A MATLAB Tutorial, Chapter 7, page 123

• Diffraction-limited applies to a system when the effect of aberrations to the pupil wavefront is smaller than the effect of diffraction at the pupil boundary

Coherent Image Simulation:

#Coherent Image Simulation
img = plt.imread('./USAF_Resolution_Chart.jpeg') # RGB image, 8-bit
img = img[:,:,0] # grayscale image, 8-bit

M, N = img.shape
L = 0.9e-3 # image length in m
du = L/M # sample interval (m)
print(f'Object is {L:.0e} m ({M} pixels), sampled every {du:.2e} m')

Nyquistf = 1/(2*du)
print(f'Thus, the highest frequency avaiable in an ideal image (Nyquist freqeuncy) is {Nyquistf:.2e} cycles/m')

# Image space co-ordinates
u = np.arange(-L/2, L/2, du) 
v = np.arange(-L/2, L/2, du)

# binary Image, type: float
Ig = np.where(img > 128, 1.0, 0.0) #img>128, replace by 1.0, otherwise replace by 0.0
ug = np.sqrt(Ig) # an ideal geometric-optics predicted image field

plt.figure(figsize=(7,7))
plt.title('Ideal image')
plt.imshow(ug, extent=[-L/2, L/2, -L/2, L/2], cmap='gray')
plt.ticklabel_format(style="sci", scilimits=(0,0))
plt.xlabel('u (m)')
plt.ylabel('v (m)')
plt.show()

print('Note that u and v are image plane co-ordinates')

Object is 9e-04 m (780 pixels), sampled every 1.15e-06 m
Thus, the highest frequency avaiable in an ideal image (Nyquist freqeuncy) is 4.33e+05 cycles/m

Note that u and v are image plane co-ordinates

lamda = 0.5e-6 # wavelength in m
wxp = 6.25e-3 # exit pupil radius in m
zxp = 125e-3 # exit pupil distance in m

fnumber = zxp/(2*wxp)
print(f'f# of the lens is {fnumber}')

f0 = wxp/(lamda*zxp) # coherent cutoff frequency
print(f'Coherent cutoff frequency with this lens is {f0:.2e} cycles/m')

if Nyquistf >= 2*f0:
    print('Imaging is in the diffraction-limited regime')
else:
    print('Not in the diffraction-limited regime: object is undersampled')

f# of the lens is 10.0
Coherent cutoff frequency with this lens is 1.00e+05 cycles/m
Imaging is in the diffraction-limited regime

#frequency coordinates
fu = np.arange(-1/(2*du), 1/(2*du), 1/L) 
fv = np.arange(-1/(2*du), 1/(2*du), 1/L) 

Fu,Fv = np.meshgrid(fu,fv)
H = np.sqrt(Fu**2 + Fv**2)/f0 < 1 # Coherent Transfer Function (Amplitude Transfer Function)

fig = plt.figure(figsize=(9,9))
fig.suptitle('Coherent Transfer Function')
ax = fig.add_subplot(projection='3d')

# Plot the surface
ax.plot_surface(Fu, Fv, H, cmap=cm.gray)
plt.ticklabel_format(style="sci", scilimits=(0,0))
fig.tight_layout()
ax.set_xlabel('fu (cycle/m)')
ax.set_ylabel('fv (cycle/m)')
plt.show()

Note that the ideal image field is blurred by the impulse response function of the imaging system, h(u,v)= FFT[H(Fx,Fy)]

i.e., Image(u,v) = Ideal_Image(u,v) convolved with h(u,v)

This is implemented in frequency-domain as follows:

Gi(Fx,Fy)= H*FFT(Ideal_Image)

So, actual image field, ui(u,v) = IFFT(Gi)

H = fftshift(H) #Coherent transfer function
Gg = fftn(fftshift(ug)) #FFT of an ideal geometric-optics predicted image field

Gi = np.multiply(Gg, H)
ui = ifftshift(ifftn(Gi))
Ii = (np.abs(ui))**2

#Simulated diffraction-limited coherent image 
plt.figure(figsize=(7,7))
plt.title('Simulated image')
plt.imshow(Ii**0.5, extent=[-L/2, L/2, -L/2, L/2], cmap='gray')
plt.ticklabel_format(style="sci", scilimits=(0,0))
plt.xlabel('u (m)')
plt.ylabel('v (m)')
plt.show()

Next, let's account for surface roughness of the object by adding a random complex exponential phase term to the field:

ug1 = np.sqrt(Ig)*np.exp(1j*2*np.pi*np.random.randn(M))

H1 = fftshift(H) #Coherent transfer function
Gg1 = fftn(fftshift(ug1))

Gi1 = np.multiply(Gg1, H1)
ui1 = ifftshift(ifftn(Gi1))
Ii1 = (np.abs(ui1))**2

#Simulated diffraction-limited coherent image 
plt.figure(figsize=(7,7))
plt.title('Simulated image of an object w/ a rough surface')
plt.imshow(Ii1**0.5, extent=[-L/2, L/2, -L/2, L/2], cmap='gray')
plt.ticklabel_format(style="sci", scilimits=(0,0))
plt.xlabel('u (m)')
plt.ylabel('v (m)')
plt.show()

The irradiance image that we get is thus impeded further by the prescence of speckles (more apparent with imshow Ii**0.5 so plotting that above)

Incoherent Image simulation:

PSF = np.abs(fftn(fftshift(H)))**2
OTF = ifftn(PSF) #left in the shifted arrangement
OTF = OTF/OTF[0,0] #normalize wrt zero-frequency value

fig = plt.figure(figsize=(9,9))
fig.suptitle('Incoherent Transfer Function')
ax = fig.add_subplot(projection='3d')

# Plot the surface
ax.plot_surface(Fu, Fv, ifftshift(np.abs(OTF)), cmap=cm.gray)
plt.ticklabel_format(style="sci", scilimits=(0,0))
fig.tight_layout()
ax.set_xlabel('fu (cycle/m)')
ax.set_ylabel('fv (cycle/m)')
plt.show()

Gg2 = fftn(fftshift(Ig)) #FFT of an ideal geometric-optics predicted image intensity
Gi2 = np.multiply(Gg2, OTF)
Ii2 = ifftshift(ifftn(Gi2))
Ii2 = np.real(Ii2) #remove residual imaginary parts
Ii2 = np.where(Ii2 <= 0, 0, Ii2) #wherever Ii2 is <= 0, replace by 0, otherwise leave as is

#Simulated diffraction-limited coherent image 
plt.figure(figsize=(7,7))
plt.title('Simulated image')
plt.imshow(Ii2, extent=[-L/2, L/2, -L/2, L/2], cmap='gray')
plt.ticklabel_format(style="sci", scilimits=(0,0))
plt.xlabel('u (m)')
plt.ylabel('v (m)')
plt.show()

Notice that the image quality in case of incoherent imaging is better than than with coherent imaging. Also, there are no ringing artifacts or speckle features in the image.